


mL = mLstar/2; 
mLub = mLstar;  % conjecture mL lies below Abel-Eberly boundary, eg, no partial replacement.

scalemL= 1/4;	% dJ_diff (below) is difference between left- and right-side J'(m) at mM. One can verify 
                % numerically that dJ_diff is declining in mL. Starting at (very) small mL, increase mL 
                % @ rate "scalemL" until dJ_diff crosses zero from above. This allows you to 
                % bracket the range in which the true mL lies.

scalemR= 1/4;   % Given an mL, we will iterate on the mR boundary: start from a low-ish mR  
                % and scale up until J'(m) crosses zero (from above). Then
                % interpolate J'() through J'()=0 to solve for mR.
mRmax = 0;      % purely a placeholder. See below.
 %
dJ_diff= 100;   
aiter  = 1;
maxaiters = 50;
while dJ_diff >0 && aiter <=maxaiters
    
    lastmL = mL;
    mL = lastmL  + scalemL*(mLub - lastmL); 


    %% Natural wastage
     % Solve for J1, J2 from value matching (VM) and smooth pasting (SP) conditions @ mL, given m_L
    muhat  = mu + (1-alpha)*s*lambda;
    rho0   = r + s*lambda;
    rho1   = rho0 - muhat;
    gamma1 = ( -(muhat-((sigma^2)/2)) - sqrt( (muhat-((sigma^2)/2))^2 + 2*(sigma^2)*rho0 ) )/(sigma^2);
    gamma2 = ( -(muhat-((sigma^2)/2)) + sqrt( (muhat-((sigma^2)/2))^2 + 2*(sigma^2)*rho0 ) )/(sigma^2);
    J2     = ((omega0/rho0) - ((1-omega1)/rho1)*((gamma1-1)/gamma1)*mL) * (gamma1/(gamma1-gamma2)) * (mL^-gamma2);
    J1     = -((mL^(1-gamma1))/gamma1) * ( ((1-omega1)/rho1) + gamma2*J2*(mL^(gamma2-1)) );
     %
     % Next, solve VM condition @ mM boundary
     % Boundary is where residm (below) crosses zero from above. (residm is 
     % difference between hiring cost c and J(m) in natural wastage region.
     % This is positive in natural wastage region but declining, reaching
     % zero at mM.) Note that residm is not necessarily monotone at m>mM.
     % Therefore, to narrow in on range where mM lies, look for minimum 
     % of residm (so long as minimum is in fact negative).
    x  = linspace(mL+.01, mL*5.0, 250)';   
    residm = c - ((1-omega1)/rho1)*x + (omega0/rho0) - J1*(x.^gamma1) - J2*(x.^gamma2);
    if min(residm) >0
        disp('ERROR: No solution to value matching problem @ mM (initialize)');
        pause;
    else
        [~,pos] = min(residm);
        if pos >1 || pos <length(x)
            % interior min
            mhat = x(pos);
        elseif pos==length(x)
            mhat = x(end);
        elseif pos==1
            disp('ERROR: No solution to value matching problem @ mM (initialize)');
        end
    end
     %
     % Note that mM must between mL and mhat --> Define x s.t. mM = mL + (mhat-mL)*(exp(x)/(1+exp(x)))
     % Note that x has support on the real line. Initial guess: mM = 1.25*mL --> 
    mM0 = 1.25*mL;
    lhs = (mM0 - mL)/(mhat-mL);
    x0  = log(lhs/(1-lhs));
    [xM, fMval, fMflag] = fzero(@(x) (c - ((1-omega1)/rho1)*(mL+(mhat-mL)*(exp(x)/(1+exp(x)))) + (omega0/rho0) - J1*((mL+(mhat-mL)*(exp(x)/(1+exp(x)))).^gamma1) - J2*((mL+(mhat-mL)*(exp(x)/(1+exp(x)))).^gamma2)), x0);
    mM = mL + (mhat-mL)*(exp(xM)/(1+exp(xM)));
    if fMflag <1
        disp('ERROR: Failed to solve for mM boundary (initialize)');
        pause;
    end 
    g  = mM/mL;
     % For later, Evaluate J'(m_M) in natural wastage region (left side derivative).
    dJ_mM_left  = ((1-omega1)/rho1) + gamma1*J1*(mM^(gamma1-1)) + gamma2*J2*(mM^(gamma2-1));



    %% Full Replacement
     % First, solve for JJ1, JJ2 using VM equs for mM and mR, J(mM) = c and J(mR) = c+C.
     % mM is already given. Solve over grid of mR. Then, solve SP equ. for mR
    Rho0 = r;
    Rho1 = r-mu;
    phi1 = ( (.5*(sigma^2)-mu) - sqrt( (.5*(sigma^2)-mu)^2 + 2*(sigma^2)*r) ) /(sigma^2);
    phi2 = ( (.5*(sigma^2)-mu) + sqrt( (.5*(sigma^2)-mu)^2 + 2*(sigma^2)*r) ) /(sigma^2);
    scriptw = ((1/(1-alpha)) - phi1) / (phi2-phi1);
     %
     % To solve for mR, we use that J'(mR) is declining as mR is pushed out further.
     % So, again, start with a low-ish mR and incrementally scale up until 
     % J'(mR) turns negative.
    quadobs = 50;
    mRobs  = 100;
    mRmax0 = 3*mM;  
    mRmax  = (aiter==1)*mRmax0 + (aiter>1)*mRmax;   % After first time through outer loop, use mRmax from last time through.
    failmR = 1;
    while failmR ==1
        aux1  = ((1-alpha)/((sigma^2)/2));
        mRvec = linspace( mM+1e-4, mRmax, mRobs)';  % vector of potential boundaries.
        mRmat = ones(quadobs,1)*mRvec';
        Gvec  = mRvec/mM;
        Gmat  = mRmat/mM; 
         %
        [nodes,wghts] = legendre(mM, mRvec, quadobs);   % quadobs by mRobs
        delta    = s*lambda ./ (1 + aux1*s*lambda*log(nodes/mM));
        integrand=  (  scriptw *((mRmat./nodes).^phi1) + ...
                    (1-scriptw)*((mRmat./nodes).^phi2) ).*(delta./nodes);
        R_mR = c*aux1 * diag(wghts'*integrand);         % replacement costs
        JJ1  = C +(c+(omega0/Rho0))*(1-Gvec.^phi2) - ((1-omega1)/Rho1)*(Gvec - Gvec.^phi2)*mM + R_mR;
        JJ1  = JJ1 ./ ( (mM^phi1)*(Gvec.^phi1 - Gvec.^phi2) );
        JJ2  = (c - ((1-omega1)/Rho1)*mM + (omega0/Rho0) - JJ1*(mM^phi1)) * (mM^-phi2);
         %
         % Now solve for mR using SP
        integrand = (phi1*   scriptw *((mRmat.^(phi1-1))./(nodes.^phi1)) + ...
                     phi2*(1-scriptw)*((mRmat.^(phi2-1))./(nodes.^phi2))).*(delta./nodes);
        dR_mR     = c*aux1*(  (delta(end,:)'./mRvec) + diag(wghts'*integrand)  ); 
        resid  = ((1-omega1)/Rho1) - dR_mR + phi1*JJ1.*(mRvec.^(phi1-1)) + phi2*JJ2.*(mRvec.^(phi2-1)); 
        if resid(end) >0
            % resid (eg, J'(mR)) is declining as mR is pushed out further. If it is
            % everywhere positive, then increase upper bound on mR and try again.
            mRmax = mRmax * (1+scalemR);
        else
            failmR = 0;
            mR     = interp1(resid,mRvec,0);
            G      = mR/mM;
            JJ1_mR = interp1(mRvec, JJ1, mR);   % recover constants given boundary, mR.
            JJ2_mR = interp1(mRvec, JJ2, mR);
        end
        
    end
     %
     % Evaluate dJ/dm @ mM in full replacement region (right side derivative):
    dR_mM = c*aux1 * (s*lambda/mM);
    dJ_mM_right = ((1-omega1)/Rho1) - dR_mM  + phi1*JJ1_mR*(mM^(phi1-1)) + phi2*JJ2_mR*(mM^(phi2-1));
    dJ_diff     = dJ_mM_left - dJ_mM_right;   
    aiter  = aiter +1;
end
 %
 %
if aiter <maxaiters
    mLlo = lastmL -.01;
    mLhi = mL + .01;        % Lower boundary lies between mLlo and mLhi.
else
    disp('ERROR: Failed to initialize lower boundary');
    pause;
end
